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Abstract 

A new cluster algorithm based on invasion percolation is described. The 
algorithm samples the critical point of a spin system without a priori knowl- 
edge of the critical temperature and provides an efficient way to determine 
the critical temperature and other observables in the critical region. The 
method is illustrated for the two- and three-dimensional Ising models. The 
algorithm equilibrates spin configurations much faster than the closely related 
Swendsen-Wang algorithm. 
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Enormous improvements in simulating systems near critical points have been achieved 
by using cluster algorithms In the present paper we describe a new cluster method 

which has the additional property of 'self-organized criticality. ' In particular, the method 
can be used to sample the critical region of various spin models without the need to fine 
tune any parameters (or know them in advance). Here, as in other cluster algorithms, 
bond clusters play a pivotal role in a Markov process where successive spin configurations 
are generated using the Fortuin-Kasteleyn || representation to identify clusters of spins for 
flipping. However, the clusters themselves are identified using invasion percolation. The new 
algorithm is closely related to the Swendsen-Wang (SW) algorithm |L| and may be adapted 
for a wide range of systems. For purposes of illustration, in this work we will consider the 
Ising model. 

Let us first recall the SW algorithm as applied to an Ising system. Starting from a spin 
configuration, satisfied bonds-those connecting spins that are of the same type-are occupied 
with probability, p((3) = 1 — where (3 = J/ksT is the coupling strength. Unsatisfied 
bonds are never occupied. Clusters of sites connected by occupied bonds are locked into 
the same spin-type and all clusters (including isolated sites) are independently flipped with 
probability 1/2. The SW algorithm samples the canonical ensemble of the spin system at 
coupling (3 and/or the random cluster (bond configuration) ensemble with parameter p. At 
T c the SW algorithm is far more efficient than single spin-flip methods because the flipped 
clusters are also critical droplets [[|. 

Here we propose using invasion percolation [|5|-|T0H to generate the bond clusters for the 



spin flips. In the usual invasion percolation, random numbers are independently assigned 
to the bonds of the lattice. Growth starts from one or more seed sites and at each step the 
clusters grow by the addition of the perimeter bond with the smallest random number. If a 
single cluster grows indefinitely on an infinite lattice, its large scale behavior is presumed to 
be that of the "incipient infinite cluster" of ordinary percolation. In particular, the fraction 
of perimeter bonds accepted into the growing cluster approaches the percolation threshold, 
p c PJTo[|. Invasion percolation is thus a self-organized critical phenomenon. 



For the present, we modify invasion percolation in two ways. First, we initiate cluster 
growth at all lattice sites. Consider this change for ordinary invasion percolation: every 
bond is initially a perimeter bond and the invasion process consists of collecting bonds in 
a given random order. Initially every site is a cluster and in most steps of the growth 
process two smaller clusters are combined into a single larger cluster. The growth process 
is terminated when some cluster first spans the system. Let / be the fraction of bonds 
accepted during the growth process and L the system size. As L — > oo, / approaches p c , the 
ordinary percolation threshold for the corresponding lattice. An intuitive argument for this 
(that can be made precise and in fact goes back to Ref. ||) is as follows: if / were smaller 
than p c the probability of a spanning cluster would vanishes as L — > oo while if / > p c , with 
probability approaching one, many disjoint subsets of the same cluster achieve spanning. 
Thus, for large systems / is close to p c . 

The second modification, which is the cornerstone of this Letter, correlates invasion 
percolation to an underlying spin configuration. As in the SW algorithm, this is done by 
allowing cluster growth along only satisfied bonds. 

The new method, which we call the invaded cluster (IC) algorithm, works as follows. 
Starting with an Ising spin configuration, S, the bonds of the lattice are given a random 
order. Correlated invasion percolation clusters are grown as described above until one of the 
clusters spans the system. (In the present implementation, a cluster is counted as spanning 
when the maximum separation in one of the d directions for some pair of points in the 
cluster is the system size, L.) After the growth process is terminated, each cluster is flipped 
with probability 1/2 yielding a new spin configuration, S' . The bonds are then randomly 
reordered and the process begins anew. 

We tested the algorithm on the nearest neighbor, ferromagnetic Ising model on square 
and simple cubic lattices with periodic boundary conditions. The computation time per 
Monte Carlo step scales linearly in the number of spins but with a somewhat larger prefactor 
than for the SW algorithm. After equilibration for 200 steps starting from the spin up state, 
statistics were collected on the energy, E and the ratio, / of the number of accepted bonds 

3 



to the number of satisfied bonds. For each system size we collected statistics for the order of 
10 4 Monte Carlo steps. We first discuss the results for two-dimensional systems with sizes 
up to 500 2 . In Fig. [I] we show the mean and median JO]] values of / plotted against 1/L. 
A linear fit through the median data extrapolates to 0.5855 compared with the exact value, 
p(/3 c ) = l-e-P c = 0.58579.... 

In Fig. 0, we plot the standard deviation of / as a function of 1/L. The solid line is a fit 
to a function of the form co + c\L~ l l 2 -\-c2L~ 1 which yields Co = —0.0014. Figures p] and § are 
thus consistent with the hypothesis that the distribution of / approaches a delta function at 
p(/3 c ) as L — ► 00 with a scaling that is given, approximately, by L -1 / 2 . The average energy 
per spin, (E)/N is shown in Fig. |3| and plotted against the inverse of the system size, 1/L. 
The best fit to the form e + t\L~ x + e 2 L~ 2 yields e = —1.706 in comparison to the exact 
result, —1.7071. The variance of the total energy divided by the number of spins is shown in 
the inset of Fig. 0. In the canonical ensemble, vax(E)/N is proportional to the specific heat 
and diverges logarithmically in L whereas here we find that this quantity diverges linearly 
in L. It is clear that the IC algorithm does not sample the canonical ensemble. 

The results for the three-dimensional Ising model up to system size 70 3 are qualitatively 
similar to the two-dimensional results except that the finite size scaling behavior for / is 
controlled by the three-dimensional Ising correlation length exponent, v ~ .63. In Fig. [I] 
we plot the mean value of / against L~ x l v . The solid line in the figure is the least squares 
linear fit to all the data which yields an intercept at .35803 in comparison to the accepted 
value fL2f , p(/3 c ) ~ .35810. The standard deviation of / vanishes with a leading behavior 
that is well fitted to L~ l l 2u . The mean energy per spin extrapolates to —2.0. 

Why does the IC algorithm work? We do not have a rigorous proof that as the system 
size goes to infinity the distribution for / peaks at p(/3 c ) or that observables such as the 
energy density converge to their limiting (infinite volume) values at criticality. Nonetheless 
we can give some heuristic arguments supporting the validity of the algorithm. 

The discussion is based on the observation that each iteration of the invaded cluster 
algorithm is identical to one iteration of the SW algorithm with p = f. Suppose we start 
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with an infinite (or huge) spin configuration that is already typical of the critical point. 
On the basis of current understanding, p((3 c ) corresponds to the threshold for the formation 
of large-scale bond clusters in the associated correlated bond percolation problem on this 
configuration. Let us now parallel the arguments given above for uncorrelated invasion 
percolation. On a finite system, in a critical spin configuration, the fraction / of accepted 
satisfied bonds will be close to p(/3 c ) when some cluster first spans the system. Hence, when 
the spin clusters get flipped, the new configuration should still be typical of criticality. It 
thus follows that if the invaded cluster algorithm starts from a critical spin configuration, 
it behaves like the SW algorithm with a temperature that fluctuates near T c and therefore 
remains near criticality. 

It is also clear that the IC algorithm moves the spin configuration toward criticality if it 
is started in either the high or low temperature phase. Suppose the spin configuration is in 
the low temperature phase. Here the portion of satisfied bonds is greater than the critical 
value and due to this relative abundance, a smaller fraction is needed to produce a spanning 
cluster than in the case of a critical spin configuration. For example, in the extreme case 
of the zero temperature configuration, one obtains / ~ p c , the ordinary percolation point 
which is of course significantly smaller than p(/3 c ). Writing / = 1 — e _/3 , this corresponds 
to an iteration of the SW algorithm at T > T c and the system is pushed toward higher 
temperature. Conversely, if the spin system is in the high temperature phase, there are not 
enough satisfied bonds and spanning will occur for / > p(/3 c ), corresponding to an iteration 
of the SW algorithm at a temperature less than the critical temperature. 

These arguments suggest that, in finite volume, the stationary distribution of the IC 
algorithm is close to (although not identical to) the canonical ensemble at f3 c and/or the 
corresponding Fortuin-Kasteleyn random cluster distribution at p(/3 c ). We will refer to the 
distribution sampled by the algorithm as the invaded cluster ensemble. Let us further 
suppose that the distribution for / becomes sharp as L — > oo and that the volume fraction 
of the spanning cluster tends to zero in this limit. It then follows that in the invaded cluster 
ensemble, the distribution functions of all local observables, e.g. spin correlation functions 
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or cluster size distributions, will converge to their infinite volume critical point distributions. 
The critical exponent, r(= 2 + /?/(/? + 7)) can be obtained from the cluster size distribution. 
Our measurements of the cluster size distribution for the two- and three-dimensional Ising 
models are consistent with the accepted values of r. Another independent exponent can 
presumably be extracted via finite-size scaling. On the other hand, we must emphasize that 
finite-volume fluctuations in the invaded cluster ensemble such as vai(E) need not have 
the same value as in the canonical ensemble and cannot be interpreted as thermodynamic 
response functions. 

We measured the normalized energy and magnetization autocorrelation functions. The 
energy autocorrelation function is defined by ((E(t) — (E))(E(0) — (E)))/vax(E) with t 
the number of iterations of the algorithm. Results for two dimensions are plotted in Fig. 
f| and compared to the SW algorithm. The energy and magnetization are almost fully 
decorrelated in a single Monte Carlo step! Results for the three-dimensional Ising model 
are similar. The negative overshoot in the energy autocorrelation is consistent with the 
negative feedback mechanism described above and the latter suggest why the algorithm is 
so fast. Consider again the example of an initial spin configuration at zero temperature: one 
iteration of the SW algorithm at (3 C yields a bond percolation configuration at p((3 c ) > p c 
which still maintains a considerable degree of low temperature order. In particular, the 
average magnetization per site is still appreciable. By contrast, the magnetization after 
one step of the IC algorithm will have essentially vanished. If the same type of reasoning 
is applied to more general configurations, the conclusion is that the IC algorithm drives a 
system to criticality faster than the SW algorithm with p = p(/3 c ). It is tempting to speculate 
that in some cases invaded cluster dynamics has no critical slowing down. 

The invaded cluster algorithm should find many uses. The extremely rapid equilibration 
time suggests it may be the best approach for high precision simulations of the critical re- 
gion of large spin systems. Using the embedding method of Ref . |2| , continuous spin models 
may be simulated. The algorithm may also be useful for first order transitions, prelimi- 
nary results for the d = 3, q = 3 Potts models indicates that the transition temperature is 
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correctly located. The IC algorithm should also prove useful for quenched random ferromag- 
netic systems where the critical temperature-which depends on the details of the disorder 
distribution-is often difficult to pin down. 

We acknowledge useful conversations with M. S. Cao during the early stages of this work. 
This work was supported by NSF Grants DMR-93-11580 and DMS-93-02023. 
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FIGURES 




FIG. 1. The fraction of accepted bonds for two and three dimensions vs L~ x l v . For two 
dimensions the mean and median values of / are shown. The solid lines are linear fits through the 
data and the known values of p((3 c ) are marked on the vertical axes with diamonds. 
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FIG. 2. The standard deviation of / vs 1/L for the two-dimensional Ising model. The solid 
line is the least squares fit to the form cq + c\L~ x l 2 + C2L~ 1 . 
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FIG. 3. (E)/N vs 1/L for the two-dimensional Ising model. The solid line is a fit to the form 
eo + e\L~ l + e2-^~ 2 and the exact infinite volume result is indicated by a diamond on the vertical 
axis. The inset shows vai(E)/N vs L. The solid line is a linear fit through the data. 
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FIG. 4. The normalized magnetization (M) and energy (E) autocorrelation functions for 
Swendsen-Wang (SW) and invaded cluster (IC) dynamics for the two-dimensional Ising model. 
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